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ABSTRACT 

We detect correlations in the cosmic far-infrared background due to the clustering of star-forming 
galaxies in observations made with the Balloon-borne Large Aperture Submillimeter Telescope, 
BLAST, at 250, 350, and 500 /im. We perform jackknife and other tests to confirm the reality of the 
signal. The measured correlations are well fit by a power law over scales of 5-25 arcminutes, with 
AI/I = 15.1 ±1.7%. We adopt a specific model for submillimeter sources in which the contribution to 
clustering comes from sources in the redshift ranges 1.3 < z < 2.2, 1.5 < z < 2.7, and 1.7 < z < 3.2, 
at 250, 350 and 500 /im, respectively. With these distributions, our measurement of the power 
spectrum, P{kg), corresponds to linear bias parameters, b = 3.8 ± 0.6,3.9 ± 0.6 and 4.4 ± 0.7, 
respectively. We further interpret the results in terms of the halo model, and find that at the 
smaller scales, the simplest halo model fails to fit our results. One way to improve the fit is to 
increase the radius at which dark matter halos are artificially truncated in the model, which is 
equivalent to having some star-forming galaxies at z > 1 located in the outskirts of groups and 
clusters. In the context of this model we find a minimum halo mass required to host a galaxy is 
log(Mmin/MQ) = 11-5+01, and we derive effective biases 6off = 2.2 ± 0.2,2.4 ± 0.2, and 2.6 ± 0.2, 
and effective masses log(Meff/M0) = 12.9 ± 0.3, 12.8 ± 0.2, and 12.7 ± 0.2 , at 250, 350 and 500 /im, 
corresponding to spatial correlation lengths of ro = 4.9,5.0, and 5.2 ± 0.7 Mpc, respectively. 
Finally, we discuss implications for clustering measurement strategies with Herschel and Planck. 
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1. INTRODUCTION 

With the discovery of the cosmic far-infrared back- 
ground (GIB, Puget et al. 1996; Fixsen et al. 1998) and 
subsequent studies from the mid-infrared to the millime- 
ter, it has been established that t he peak epoch of star 
formation lies between 1 ^ z ^ 3 (jDickinson et al.ll2003l : 

[Hopkins 2004). 

Where star formation occurs with respect to the under- 
lying dark matter distribution is less well understood. In 
the local Universe, the sites of most active star-formation 
occur far from the den sest environments , a property 
which reverses at z ~ 1 (l Elbaz et alll2007 r). This points 
to the environment contributing to the mechanisms that 
trigger or quench star formation. 

Measurements of the clustering of star-forming galax- 
ies on large and small scales can be used to directly relate 
star formation to the environment. In the regime of lin- 
ear growth of structure, the clustering amplitude relative 
to that of the underlying dark matter is described by a 
bias parameter, 6, which describes how strongly star for- 
mation traces the underlying dark matter. On smaller 
scales, the distribution of galaxies hosted by a dark mat- 
ter halo is described by the halo occupation distribution 
(jPeacock fc Smit3l2000D . It contains information on the 
abundance of star-forming sources within individual ha- 
los. Targeting the submillimeter band is particularly ef- 
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ficient for observing star formation directly. The sub- 
millimeter background results from the thermal emission 
of interstellar dust in high-redshift star-forming galax- 
ies, which is heated by optical and ultraviolet radiation 
from stars and to a le sser extent active galactic nuclei 
(see iBlain et aD l2002| ) . A substantial effo rt has been 



(IScott fc WhitelllQQaiHaiman fc Knox l[2000l: IKnox eFall 
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devoted to surveys of these galaxies (e.g., ISmail et al 



IQQTtlHughes et al.lll998l:lBarger et aLlll998i:lBorvs et al 



20031) : thus, in principle, a clustering signal could be mea- 
sured from these sources directly. 

However, direct measurement of the clustering proper- 
ties of resolved submillimeter galaxies has been elusive. 
Due to limited mapping speeds, the areal coverage of 
even the most ambitious submillimeter surveys has been 
relatively small (e.g., SHADES mapped ap proximately a 
quarter square degree at 850 /im; iCoppin et al. 2006) . In 
addition, due to steeply falling counts and modest resolu- 
tions of single-dish submillimeter telescopes, source con- 
fusion has made it difficult to resolve any sources other 
than those with very high signal-to-noise. These sources 
span a relatively wide redshift range, rough ly 1 < z < 4, 
peaking at z ~ 2.4 (jChapman et al.l 120051 ). which has 
the effect of washing out the angular clustering sig- 
nal, further compl icating measureme nts. This difficulty 
was confirmed by IScott et al.l (|2006f) . who re-analyzed 
all the SCUBA fields and found tenative evidence of 
strong angular clustering, but with errors too large to 
adequately constrain the spatial correlation length. Fur- 
thermore, relatively large beams have ma de it difficult 
to efficiently ideiitify d irect counterparts (jBarger et al.l 
119991 : ITvison et al]|2000D ii i order to ob t ain sp ectroscopic 
or photometric redshifts. IBlain et al.l (|2004f ) attempted 
to measure the spatial clustering properties combining 
73 sources with spectroscopic information, but again 
were only able to tentatively measure the clustering 
length. Other SCUBA galaxy clustering measurements, 
so me tentativel y detecting clustering, h a ve be en made 
by IWebb et all (120 03) and IBlake et all ((200l) : mean- 
while I Almaini et aLI (j2003) claim to have found evidence 
for strong angular clustering bet ween X-ray and su bmil- 
limeter populations, although Borvs et al.l (|2004f l were 
unable to confirm the result using another, seemingly 
less biased, estimat or. Using a nearest-neighbor analysis, 
iGreve et all ()2004D find that most significant MAMBO 
(1.2 mm) sources come in pairs, separated by ~ 23 arcsec. 
Put together, these limitations have made it extremely 
difficult to measure the clustering signal robustly. 

To study the clustering properties of submillimeter 
galaxies it is more powerful to consider the statis- 
tics of the unresolved CIB, which contains the full 
intensity, rather than working from a limited cata- 
log. In other words, instead of measuring a correlation 
among numbers of galaxies, use the background fluc- 
tuations of the total intensity to measure c orrelations 
among brightn e sses o f galaxies. iDevlin et al.l (2009) and 
iMarsden et al.l (|2009f ) have demonstrated that the CIB 
is composed of emission by discreet sources. Since sub- 
millimeter galaxies are optically thin, this signal will be 
proportional to the total star formation rates of those 
sources. Correlations in the CIB will have a contribu- 
tion in excess of white noise - which arises from Poisson 
sampling of a background made up of a finite number 
of sources - in the presence of clustering, with an am- 
plitude that should be detectable with current surveys 



_ , ■ .Magliocchetti et a.1.1 12 001: P errotta et al.l 120031: 
Amblard fc Cooravl l2007nNegrello et al.ll2007D. Initial 
attempts to detect correlations bv lPeacock etaH (120001. 
for the Hubble Deep Field observed by SCUBA at 
850 /im, and bv lLagache fc Pugetl (pOOOl) for a 0.25 deg^ 
ISO field at 170 /xm, were only able to measure a signal 
consistent with th e Poiss on co ntribution. Mo r e rece ntly, 
iGrossan fc SmootI ([20071 ) and iLagache et"all (|2007D re- 
ported the weak detection of a clustering component in 
160 /im data from ^ 9 deg'^ Spitzer fields. 

In this paper we report the detection of correlations 
in the submillimeter part of the CIB due to the cluster- 
ing of star-forming galaxies, in a 6 deg^ field centered 
on the Great Observatories Origins Deep S urvey South 
field (GOODS-South; Giava lisco et al.|[200l . These data 
were collected by the Balloon- borne Large Aperture Sub- 
millimeter Telescope fBLAST: [Devlin et a l. 2009), which 
is designed to bracket the peak of redshiftcd thermal 
emission from dust by observing at 250, 350, and 500 /im. 
Operating above most of the atmosphere, BLAST is able 
to make observations in bands which are difficult or im- 
possible to observe from the ground. A detailed descrip- 
tion of the instrum ent a nd calibration can be found in 
iPascale et all ([200l and lTruch et~aI1 ((20091 ). 

This paper is organized as follows. In Part I we de- 
scribe how we make the measurement - from map prepa- 
ration to power spectrum calculation. We address each 
contribution to the total power spectrum, and how they 
are removed to uncover the clustering signal. We show 
that the observed spectra are consistent across BLAST 
bands and correspond to a modest bias. In Part II we 
assess the plausibility of the detection by fitting models 
to the data: beginning with a simple linea r bias model, 
followed by a more detailed halo model ([Mo fc Whitd 
119961 : iCoorav fc She th 2002 ) and halo occupation dis- 
trihution ([Peacock fc Smith l2000() . When required we 
adopt the concordance model, a flat ACDM cosmology 
with = 0.2 74, n^ = 0.726 gp = 70.5 km s"! Mpc-\ 
and cTg = 0.81 ([Hinshaw et al.ll2009[ ). 

PART I: MEASURING CORRELATIONS IN THE CIB 

2. BACKGROUND CORRELATIONS AND THE 
CORRELATION FUNCTION; OVERVIEW 

Galaxy clustering can be expressed in a number of 
ways, the most common being the two-point correlation 
function, w{9), which measures the number of pairs at 
a given distance in excess of what would be expected of 
a Poisson distribution. Alternatively, the clustering of 
galaxies can be expressed as a power spectrum in excess 
of Poisson noise, P{ke), which can be expressed in dimen- 
sionless units as the fractional variance per logarithmic 
increment of wavenumber fsee lPeacocklll999l ). i.e.. 



A?, = 



(1) 



where kg is the angular wavenumber, which is also known 
as a in the literature, and is expressed in inverse angular 
scale as kg = 1/A. It is related to the multipole index, 
i,hy i ~ 2TTke. The tildes over and P{ke) denote 
that these quantities refer to galaxy locations rather than 
intensities. 

Naturally, the correlation function and the power 
spectrum of galaxy clustering are related; they form a 
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Hankel transform pair. Explicitly, for small surveys, 
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For small separation correlations of local galaxies, angu- 
lar clustering is often described as a power law, w{9) = 
(9/9o)~'^, were e ~ 0.8 is the canonical slope (e.g., 
iGiavalisco et al.lll998D. The p ower spectrum analog for 
small areas (see lPeacock et al.ll2000( ) is 



Alike) = {2TTkeeor2^~ 



. r(l-e/2) 

r(6/2) ' 



(3) 



which is equal to 3.35(fce6lo)"-« for e ~ 0.8. 

As previously stated, this holds for correlations of 
galaxy locations. On the other hand, the power spec- 
trum of background fluctuations, which is what we cal- 
culate, measures correlations of galaxy intensities. The 
redshift distribution of the cumulative flux contributed 
by the background sources is represented as 



dS 

dz 



dN 

S*— — dS, 
db dz 



(4) 



where dN/ (dSdz) is the number density of sources per 
unit flux density and redshift interval. The measured 
power spectrum of background fluctuations, P{kg), is 
the 2-dimensional, flux weighted, projection of the 3- 
dimensional galaxy clustering spectrum, Psuike). For 
kg ^ 1 and a flat cosmology, their relationship can be 
approximated as 



P{ke) 



P3D(2.fc./.(z),z)(f(z)) 



dz, 
(5) 

(e.g.. iTegmark et aLll2002l) . where x{z) is the comoving 
radial distance and dVc{z) is the comoving volume ele- 
ment, i.e., dVc{z) = x{z)'^dx/dz (assuming a flat cosmol- 
ogy), where dN/{dSdz) is the number density of sources 
per unit flux density and redshift interval. 
P{kg) can expressed in dimensionless units as 



^ 2TTklP{ke)/ll, 



(6) 



where Ii, is the intensity of the background in Jy. Al- 
though A^^ and A^^ have the same units and similar 
meaning, they differ because Equation [6] deals with an 
intensity weighted power spectrum. 

In addition to the clustering signal, the total power 
spectrum has contributions from instrumental noise, 
Poisson noise from individual background galaxies, and 
cirrus emission. We will address each contribution indi- 
vidually. 

3. METHODS 
3.1. Map Preparation 

BLAST observed a wide 8.7 deg'^ patch, centered 
on the GOODS-South field (3h32"35", -28°15'; here- 
after BGS-Wide), with mean 1-a sensitivities of 36, 
31, and 20mJybeam~ at 250, 350, and 500 /im, as 
well as a deep, nested field of 0.8 deg^, centered on 



(3'^32"^30", -27°48'; hereafter BGS-Deep) with mean 1- 
(T sensitivities of 11, 9 and 6 mJy beam""'^, respectively^. 
A 6 deg^ region, centered on BGS-Wide, is selected from 
the map, because of its uniformity in observed depth. 
The BLAST bolometers are prone to drifts on timescales 
greater than ~ 10 seconds. To retain as much large an- 
gular scale signal as possible, a fast scan rate is preferred 
because larger scales will survive the high-pass filtering, 
designed to remove noise below the 0.1 Hz 1// knee. 
For this analysis we use only the data from the wide 
region of the map, where the scan rate is O.ldegs^^, 
and the r.m.s. of the maps is 2.4, 1.9, l.OMJysr"^, 
at 250, 350, and 500 /zm, which is applicable to calcu- 
lating uncertainties of point-source flux densities. The 
data for the nested deep region, whose scan rate is only 
0.05degs~^, are not included. Large-scale noise is re- 
moved by high-pass filtering the time-streams at 0.2 Hz. 
Correlated noise - a drift of multiple detectors in unison - 
is not removed because it would inevitably suppress large 
scale signal as well. Instead, a cross-correlation of a sub- 
set of the maps is used to remove large-scale correlated 
noise (described below). Although fully optimized map- 
make rs are available (e.g., SANEPIC; ,Patanchon et al.l 
|2008| ). the long time to convergence (typically 24 hours 
on 10 processors for this particular map) makes it im- 
practical to use them for our Monte-Carlo simulations. 
Full analysis with SANEPIC maps, including the nested 
deep field, will be the subject of a future p aper. Here we 
instead use OPTBIN (jPascale et al.l 12003 ) a fast, naive, 
map-maker whose transfer function is calculated using a 
Monte-Carlo simulation (see § 13.21 and Figure [T] for de- 
tails), and we use SANEPIC maps for consistency checks. 
Parallel power spectrum analysis with both map-makers 
show agreement well within the simulated uncertainties. 

3.2. Power Spectrum Calculation 



The map intensity, 5map, can be written as 

^map = (T ® [S-sky ® B + N]) W, 



(7) 



where we use ® to represent a convolution, 5'sky is the 
true sky surface brightness, T is the transfer function of 
the map-maker, B is the measured instrumental beam, 
N is the instrumental noise, and W is the 'aperture 
function', which is zero beyond the region of interest. 
The autocorrelation of a map will contain a contribution 
from detector noise. To suppress this instrumental noise, 
cross-power spectra are taken among a set of four maps 
which are made by dividing the time-stream into four 
roughly equal parts and then making four separate maps 
(hereafter referred to as sub- maps). The timestreams, 
which are made up of numerous chunks, are divided into 
every fourth chunk (e.g., 1, 5, 9, ... , and 2, 6, 10, . . . , 
etc.), so that the sub-maps have as similar coverage as 
possible. The number of sub-maps chosen maximizes the 
number that can be made while maintaining uniformity 
in hits, retaining some cross-linking, and avoiding holes 
in the maps. The r.m.s. of the resulting sub-maps is 4.6, 
3.6, and 2.0MJysr-i, at 250, 350, and 500 /xm. In the 
cross-spectrum, noise which is uncorrelated between sub- 
maps averages to zero. Consequently, the spectrum does 

^ BLAST maps and catalogs are publicly available at 
http: / /www. blastexperiment. info 
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Fig. 1. — Transfer functions calculated with a Monte-Carlo sim- 
ulation involving 500 mock-maps observed with the BLAST simu- 
lator. The lines are identical beyond the scale of the beams, which 
implies that the map-maker is linear, as we would expect. 

not depend on modeling the potentially complicated or 
non-stationary noise. 

We prepare the maps before calculating the power 
spectrum by remo ving their m eans, apodizing them with 
a Welch window ('Presd [2b02l . chapter 13.4), and zero- 
padding them with a width on each side equal to half the 
map. The cross-correlation two-dimensional power spec- 
trum of each pair of maps is calculated. The azimuthal 
average of the amplitudes (which in two-dimensional k- 
space appears to be isotropic), is taken to find the one- 
dimensional power spectrum, P{k). The resulting spec- 
tra are averaged and divided by the power spectrum of 
the beam and the transfer function. The transfer func- 
tion is calculated with a Monte-Carlo simulation from 
simulated maps made with the BLAST simulator^, and 
is shown in Figure [1] At angular scales which are large 
compared to the high-pass filter (< 0.03 arcmin"^) and 
approaching the cut-off of the beam (> 0.9 arcmin"^), 
the transfer function is unreliable. We are interested in 
the scales bracketed by these limits. 

3.2.1. Jackknife tests 

We have performed jackknife tests in which a cor- 
relation is calculated between two distinct difi'erence 
maps. Specifically, we take the cross-correlation of: (sub- 
map 1 — sub-map 2) and (sub-map 3 — sub-map 4). If the 
cross-correlation measures only signal, then taking a dif- 
ference between sub-maps should cancel sky signal and 
result in a cross-correlation power spectrum consistent 
with zero. And indeed, our results are consistent with 
zero, in the range of interest, in all three bands. 

As a further check, the cross-power spectrum is com- 
pared to the difference of the auto-power spectrum (of 
a map made from the entire timestream) and the mea- 
sured noise. The noise is estimated from odd-even pixel 
jackknife maps, and is approximately white over the an- 
gular scales of interest. The resulting two spectra are in 
excellent agreement. 



^ We simulate the BLAST pipeline by 'observing' an input map, 
creating a set of time-streams, filtering them, and making a map 
from the filtered timestreams. 



3.2.2. Poisson Noise 

Poisson (or shot) noise arises from the finite number of 
galaxies per unit area. It can be calculated analytically 
from the source counts as 

Psi^ot ^ J^°° S^^dS. (8) 

Alternatively, Poisson noise levels can be esti- 
mated from Monte-Carlo simulations using mock-maps 
which are populated with uncorrelated sources whose 
fluxes are drawn from the measured BLAST counts 
(jPatanchon et al.] [2009). This has the added advan- 
tage that realistic uncertainties can be calculated as well. 
Since the counts are so steep, care must be taken to re- 
produce the bright end of the counts faithfully, so that 
extremely rare bright sources do not unrealistically ap- 
pear in the realizations. We find that the two methods 
of estimating the level of shot noise agree within the 1-cr 
uncertainties. 

Due to the steep nature of the source counts, the 
Poisson noise is dominated by th e contributions of the 
fainter population. Furthermore, IMarsden et al.l (j2009f ) 
find only 15% of the total sky intensity is associated with 
a 3-a catalog. We find that removal of only the bright- 
est sources results in an approximately 5% reduction in 
Poisson noise at 250 /Ltm, and that more aggressive cuts 
lead to removal of correlations beyond just Poisson noise. 
We find behavior consistent with this in our simulations. 
Therefore, we subtract only 5 sources above 500 mJy 
at 250 /im, 2 sources above 400 mJy at 350 /im, and no 
sources at 500 /.tm. To subtract sources, first we make a 
source list by performing a noise- weighted convolution of 
the maps with the effective BLAST point-spread func- 
tion (PSF) and identify local maxima in the smoothed 
map. We then subtract a scaled effective PSF with am- 
plitude taken from the source list. We perform the same 
operation on our mock-maps, from which we calculate 
Poisson levels of 11.4 ± 1.0 x 10^,6.3 ± 0.5 x 10^, and 
2.7 ±0.2 X lO^Jy^sr-i at 250, 350, and 500 /im. These 
are shown as dashed lines in Figure [2l 

3.2.3. Estimates of Galactic Foregrounds 

iG autier et all ()1992f) show that the power spectrum of 
Galactic cirrus can be approximated by a power law, 

Pcirrus(fce) = Po (j^^ , (9) 

where kg is the angular wavenumber in inverse ar- 
cminutes, and Pq is the power spectrum value at 
ko = 0.01 arcmin^^. We measure the cirrus com- 
ponent at 100 /im from the cross-correlation of the 
three co- added IRIS (HCON^) map s (reprocessed IRAS: 
IMiville-Deschenes fc Lagachell200"5l) for a ^ 15 deg^ re- 
gion surrounding the BGS-Wide field. The amplitude 
of the observed power spectrum has a contribution from 
cirrus emission which is highly variable on the sky. The 
GOODS region was specifically chosen because it is 
a low cirrus region, with a mean intensity of 1.39 ± 
0.18 MJy sr-i. At scales 0.008 < kg < 0.03arcmin~\ 

^ HCON refers to each survey. For more information and maps 
see http: / / www.cita.utoronto.ca/ ~mamd/IRIS/IrisOverview.html ^ 
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Fig. 2. — Power spectra of 6 deg^ regions selected from the BGS- 
Wide maps at 250, 350, and 500 /im are shown with 1-cr uncer- 
tainties. Color-corrected fits to galactic cirrus measured from IRIS 
100 fim maps are shown as dotted lines with shaded region repre- 
senting 1-cr uncertainties, sloping down from the left side (500 fim 
cirrus is too low to make it onto the region plotted). Poisson 
noise contributions are found with Monte-Carlo simulations, and 
are shown as horizontal dashed lines with similar error regions. 
Scale invariant, k~'^, power spectra are shown as dot-dashed lines. 
Vertical dashed line at kg = 0.33 arcmin"^ represe nts the scale at 
which the variance in the FIDEL catalog — which iMarsden et al.l 
H20091 show resolves most of the CIB — is no longer Poisson noise 
dominated. For angular scales greater than 0.33 arcmin"^, signal 
in excess of Poisson noise is attributed to the clustering of star- 
forming galaxies. 

which is larger than the scales probed by BLAST, we 
find that the cirrus is well approximated by a power law, 
Po = (0.47 ± 0.18) X 106 Jy^sr-i and a = 2.91 ± 0.11, 
values which are similar to those found in low c irrus 
emission regio ns measur e d bvjLagache et al.l ()2007D and 
[Miville-Deschenes et al.l ()2007l) ^ We assume the power 
spectrum continues to smaller angular scales, and scale 
it to the BLAST bands using the average dust emission 
color (/blast/^ioo)^j which is found by assuming cir- 
rus emission behaves as a modified blackbody, v^B{v)^ 
where B(v) is the Planck fu nction, and (3 is the emissivity 
index (jPraine fc Led [198^ . The scaled power spectrum 
and errors are calculated with a Monte-Carlo simulation, 
varying temperature ( 17.5 ± 1.5 K) and (3 (1.9 ± 0.2) 
(|Boulanger et al.lll99 6D. The resulting power law approx- 
imations and uncertainties are illustrated as dotted lines 
in Figure [H which show that we are negligibly affected 
by Galactic cirrus on all recovered angular scales. 



3.2.4. Uncertainties 

To estimate uncertainties in our angular power spectra, 
clustered signal plus noise simulated maps are analyzed 
with the same pipeline as the astronomical data. We 
include clustering in the simulations in case the ampli- 
tude or shape of the transfer function depends on the 
inp ut, but we fou i id tha t not to be the case. We fol- 
low |AlmainretI2] (|2005l . see appendix for algorithm) to 
introduce correlations to the simulated maps, with an 
angular correlation length, = 3.0". This angle was 
chosen from the measured upper limit of iPeacock et al.l 
(|2000f ). Realizations with stronger clustering clearly do 
not look like BLAST maps. 

4. BASIC CLUSTERING RESULTS 

Figure [2] shows an unambiguous signal in excess of 
Poisson noise on scales of 0.04 — 0.2 arcmin"^ which can- 
not be explained by Galactic cirrus, and which we inter- 
pret as correlations from clustered star-forming galax- 
ies. The vertical dashed line indicates the angle at 
whi ch the distribut ion of 24 /im selected FIDEL galax- 
ies (jMagnelli et al.r 2009) begin to show variance in ex- 
cess of Poisson. Ind eed, the similarity to Figure 3 of 
IMarsden et"!!! ()2009f ) is striking. 

The CIB has an amplitude I^^^ measured to be 0.71 ± 
0.17, 0.59 ± 0.14, and 0. 38 ± 0.10 MJysr-^ at 250, 350, 
and 500 ^m, respectively ([Marsden et al.|[2009f) . Figure[3] 
shows the clustering component of the power spectrum 
normalized by /^^^ in the three BLAST bands, where 
the 250 /im and 500 /im data have been displaced slightly 
for visual clarity. The contributions of cirrus and Pois- 
son noise have been subtracted, and the data rebinned in 
logarithmic intervals. These results are also listed in Ta- 
ble [1] We use the BLAST estimates for the CIB because 
they are the most precise estimate available of the CIB 
in these wavelength bands. Doing so has the additional 
benefit in this case that calibration uncertainties com- 
pletely vanish in the ratio. The fit to a single power law, 
and also the agreement in amplitude across the BLAST 
bands, once the power spectra are normalized to the sky 
intensity, are excellent. 

The relative variance of the CIB, , formed by divid- 
ing P{kg) by 27rfcg, is shown in the bottom panel of Fig- 
ure [3l where the 250 /zm and 500 /xm data have been dis- 
placed slightly for visual clarity. (Amplitudes have been 
additionally converted from arcmin"^ to steradian"^ in 
the figure and in the discussion below.) The best fits of 
Equation [3]) are shown in Table [2l The data are consis- 
tent (to within 1—a) with Aj. = Constant , with the same 
amplitude in all three bands , even taking into account 
that the error bars shown contain a large component of 
uncertainty which is common mode between channels, 
and thus overestimate the anticipated scatter. The best- 
fit amplitude for a power law with slope of —2 is shown 
as a dotted line in the bottom panel, corresponding to 
Akf, = SI /I = 15.1 ± 1.7%. This is directly analogous to 
the square root of the 'band-power' measured in Cosmic 
Microwave Background (CMB) anisotropy experiments, 
where typically ST/T ^ 10^^. From this point of view 
the CIB is much clumpier than the CMB, as one would 
expect since the galaxies are observed after a much longer 
period of linear growth. The dotted line in the upper 
panel of Figure [3] is 27rfcgA? , corresponding exactly to 
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kg 

(arcmin" 



BLAST 250 

(steradian^-'-) 



0.044 
0.070 
0.113 
0.183 



(1.39 ±0.52) X 1(F^ 
(5.94 ± 2.02) X 10"* 
(2.99 ±0.95) X 10"* 
(9.81 ± 7.51) X lO-'^ 



BLAST 350 
(steradian"-'-) 



(1.33 ±0.47) X 10^ 
(4.98 ± 1.83) X 10"* 
(2.62 ± 0.90) X 10"* 
(9.94 ± 6.92) X IQ-^ 



BLAST 500 
(stcradian^^) 



(2.47 ±1.17) X 10^ 
(5.71 ±2.06) X 10-8 
(2.98 ± 1.00) X 10-8 
(6.91 ±4.62) X lO"'^ 



CIB NORMALIZED CLUSTERING POWER SPECTRA, P{kg)/l'^ 



TABLE 1 

The ERRORS DO NOT INCLUDE UNCERTAINTIES IN THE CIB. CIB VALUES ARE 
LISTED IN §[4] 



BAND 


00 (arcsec) 


e 




250 /im 


0.017 ±0.020 


0.27 ±0.19 




350 fim 


0.008 ± 0.006 


0.26 ±0.19 




500 Atm 


0.0 


0.0 






TABLE 2 




CO 



Best-fit values obtained for 9o and e. 

the fit in the lower panel. 

The values of the parameters e and (see Equation[3]) 
which best fit the data are shown in Table [21 The large 
uncertainty quoted for 9q is due partly to the awkward- 
ness of this parameterization near to zero slope. The 
bottom panel of Figure [3] shows the best-fit for 250 /im 
data as a dashed line. Such a small despite signifi- 
cant power in excess of Poisson noise, requires that the 
sources which make up the CIB are distributed over a 
wide range of redshifts, and has implications for future 
clustering measurement strategies (e.g., with Herschel or 
Planck), which we discuss in §[TTJ 

We have calculated the correlations between different 
bands (e.g., 250 x 350) using the same pipeline and set 
of sub-maps. The cross-spectra are normalized by the 
square root of the auto-spectra of the two bands, so that 
the final curve would be unity at all scales for identical 
maps, and zero at all scales for two completely different 
maps. Results show that the cross-correlations are 0.95 ± 
0.06, 1.06±0.09, and 0.92±0.04, for 250 x 350, 350 x 500, 
and 250 x 500, respectively, over the range of angular 
scales 0.04 < kg < 0.5 arcmin"^. Neighboring bands 
are more correlated with each other than are the 250 
and 500 /im bands. While we find the same spectrum in 
all three bands, the phases are different, as one would 
expect if the three BLAST bands have different selection 
functions and sample the galaxies in the CIB at different 
redshifts. While the cross-band correlation provides a 
powerful tool for testing the redshift distributions and 
spectral energy densities of source population models, 
a more detailed analysis is required before making any 
strong conclusions. This will be the subject of a future 
study using the SANEPIC maps and including the BGS- 
Deep data. 

We have measured the variance projected onto two di- 
mensions, but galaxies are of course distributed in three 
dimensions. Knowledge of the redshift distribution of the 
galaxies in the CIB allows interpretation of these power 
spectra in terms of a bias factor quantifying the compar- 
ison to cold dark matter (CDM) spectra, which we find 
in the next section. As the angular scales we probe get 
smaller, non-linear growth eventually sets in, which we 
examine by fitting halo models to these spectra later in 
the paper. 




0.03 



0.05 



0.10 
kg (arcmin"') 



0.25 



Fig. 3. — Top: Power Spectra of the clustering component, after 
removal of cirrus and Poisson noise, and normalized by {l'^^^)^ ■ 
The best-fit power spectrum, proportional to , is shown as a 
dotted line. Bottom: A| shown along with best-fit power spectra 

(horizontal dotted line); as well as best-fit parameters to Equa- 
tion |3] (dashed line) for 250 /xm data. The 250 and 500 /^m points 
are offset horizontally for clarity, by factors of -0.025 and +0.025, 
respectively. Clearly, the power spectrum signal of clustered star- 
forming galaxies is well fit by a power law power spectrum propor- 
tional to k~^ . 

PART II: MODEL FITTING 
5. INTRODUCTION 

To interpret our detection of correlations requires com- 
parison to an underlying model whose parameters the 
data constrain. Such a model could contain details of the 
complete source population, including number counts, 
i.e. intensity distributions, and redshift distributions, as 
well as a framework for describing the linear and non- 
linear clustering regimes. This latter part might involve, 
for example, a halo occupation distribution which ac- 
counts for the galaxy distribution in a given dark matter 
halo as a function of lumino sity (i.e., the so-called con- 
ditional luminosity function; lCooravl[2006D . For a back- 
ground of primarily unresolved sources, the conditional 
luminosity function model is an improvement on the sim- 
ple halo occupation distribution; however, its increase in 
complexity comes with an increase of free parameters. 

Given the very good fit to a constant A^^ , we first ex- 
plore the physical meaning of the BLAST power spectra 
in the context of a purely linear model. We assume that 
the galaxies which comprise the CIB have a power spec- 
trum which is scaled from the power spectrum of dark 
matter at every redshift, P{k) = b'^PBuik)- Because 
PuM is redshift dependent, estimating requires knowl- 
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Fig. 4. — Redshift distribution of the cumulative flux con- 
tributed by th e background s ources at the BLAST bands, ac- 
cording to the ILagache et al.l 1)20041 ) model. The dashed (blue) 
and dot-dashed (red) curves are for 'regular' and 'star-forming' 
IRAS galaxies respectively, while the solid line is the total. 
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Fig. 5. — Contribution to the angular power spectrum from dif- 
ferent redshift slices, inferred from Figure |4] at specific angular 
scales probed by the data. The different curves correspond to the 
following scales: dotted (green) - kg = 0.04 arcmin"^; solid (red) - 
kg =0.1 arcmin"^; dashed (blue) - kg = 0.18 arcmin"^. Whereas 
the redshift distribution of sources (Figure|4)l has a significant con- 
tribution from z < I, the contribution from those sources to the 
angular power spectrum is negligible. 

edge of the redshift distribution of the galaxies which 
comprise the BLAST sign als. To find this distribution 
we adopt the model of La gache et al.l (|2004D . described 
in Section [HI It is worth noting that since we are really 
measuring an emission-weighted bias at rest-frame far- 
IR wavelengths, there will be some degeneracy between 
the value of b and changes in the redshift distribution or 
far-IR spectral shapes assumed for the sources. 

We extend the fit into the non-linear regime to study 
whether the BLAST data can constrain parameters of 
the halo occupation distribution. We can check whether 
our correlations are consistent with the model (e.g., by 
judging whether the same bias fits all three BLAST wave- 
bands, and that the cross-band correlations of simula- 
tions made with the model agree with the data), but we 
will not explore how the model might be improved. It is 
important to understand that if the distribution of source 
redshifts or counts were different, then we would infer 
a different bias level. However, since the model we've 
adopted agrees with a large body of multi-wavelength 
observations, we are confident that these prescriptions 
are a reasonable first attempt. 

6. SOURCE POPULATION MODEL 

Knowledge of the full redshift distribution of BLAST 
sources would allow us to estimate the redshift distribu- 
tion of the measured clustering signal, i.e., the redshift 
range probed by the power spectrum of correlations due 
to clustering descr ibed by the P x dS ldz distribution. 
iDevlin et al.l (|2009f) . and in more detail iPatanchon et al.l 
l200g), fi nd the n umber counts of the BLAST sources, 
and .Pascale et al.l ((2b09j ) divides the redshift distribution 



roughly into four redshift bins. Precisely understanding 
the redshift distribution of those sources is a work in 
progress (e.g., Chapin et al., in prep.) . In th e meantime 
we adopt the model of ILagache et al.l ()2004[ ) - a model 
which approximates the counts and redshift distributions 
of the populations expected to make up the CIB - to de- 
scribe the underlying source population^. Specifically, it 
is a phenomenological model which extrapolates the local 
60 /im luminosity function of IRAS sources (divided into 
'regular' and 'star-forming' components) to longer wave- 
lengths, assuming a set of spectral energy distribution 
templates. The extrapolation to higher redshift is param- 
eterized as a mixture of luminosity and density evolution, 
and is constrained to be consistent with most of the avail- 
able data on source counts, redshift distributions and the 
far-infrared background intensity. It becomes less reli- 
able with greater extrapolation to longer wavelengths, 
and therefore should be considered a n approximation; 
howe ver, it does reproduce the BLAST (jPatanchon et al.l 
[2001 counts at a level which is sufficient for our pur- 
poses. Furthermore, since the clustering signal is dom- 
inated by the faint source population, it is critical that 
the model is consistent with the level of the background 
at the BLAST wavelengths, and that the cross-band cor- 
relations of simulated maps agree with those measured, 
which we find to be the case. 

Most of the clustering signal is coming from relatively 
high redshifts. The medians of the redshift distributions 
in Figure [5] are z = 1.61, 1.88 and 2.42, at 250, 350 and 
500 /im, respectively, and the upper and lower quartiles 
of the distributions are z =(1.3, 2.2), (1.5, 2.7), (1.7, 3.2), 
respectively. At the representative scale of 0.1 arcmin^^ 
(red solid line in Figure O, 95% of the background orig- 
inates from sources 1.1 < z, 1.2 < z, and 1.4 < z in the 
three bands. 

7. LINEAR BIAS MODEL 

We first carry out a simple fit to a scaled version of the 
linear theory power spectrum of the dark matter Pdm- 
As can be seen in Figure [6l a simple biasing prescription 
provides a good fit to the BLAST data. The required bias 
levels are 6 = 3.8 ± 0.6, 3.9 ± 0.6, and 4.4 ± 0.7 at 250, 
350 and 500 fim, respectively, with a reduced Xmin 0-4 
(with 10 degrees of freedom) in all three bands. 

More detailed modeling could be attempted. In prin- 
ciple the cross-band measurements and wavelength de- 
pendence of the measured correlation amplitudes could 
be used to es timate the v a riatio n of bias with redshift, as 
discussed bv IKnox et al.l (|2001[ ). However, we leave this 
to a future study. 

8. HALO MODEL 

As illustrated in FigurelHl ^ simple biasing prescription 
provides a good fit to our data. The main drawback of 
this simple fit is that it may not realistically account for 
the 1-halo, nonlinear clustering component. Our data 
are expected to bracket the physical scales correspond- 
ing to the transition from linear to non-linear clustering 
regimes, and though these two contributions appear to 
have combined to look very much like a scaled linear 

D ocumentation and IDL files can be found at 
|http: / /www.ias.u-psud.fr/irgalaxies/model.php^ 
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Fig. 6. — Power spectrum of background correlations (circles with error bars) overlaid with the best-fit linear bias (solid line), as well as 
predictions obtained for halo models different values of rcut- See Table [3] for the fit parameters. Although the power spectra shown here 
have different amplitudes in the three channels, P(kg)/Iy is the same in all three bands, as shown in Figure[3] The data are fit best by a 
model which includes a linear term only. 



bias, to conclude that there is no contribution from a 1- 
halo term may be unphysical. Therefore, in this section 
we use a particular implementation of the 'halo model', 
which assigns galaxies to halos as a function of mass, and 
consequently probes into the territory of non-linear fluc- 
tuations, to explore how the 1-halo term could contribute 
power on small scales. This also allows us to discuss our 
results in the context of other measurements of galaxy 
clustering. 

The halo model of large scale structure has proven to 
be a powerful tool for describing the c lustering proper- 
ties o f cosmic objects (for a review, see iCoorav fc ShethI 
[200l) . Its main ingredient is the p arameterization of 
the h alo occupation distribution (HOD. Peaco ck fc SmithI 
120001 ) ■ which describes how galaxies populate dark mat- 
ter halos as a function of halo mass. The power spectrum 
of galaxies is written as the sum of two components: the 
1-halo term, Pih, which describes pairs of objects within 
the same dark matter halo, and the 2-halo term, P2h, 
which accounts for pairs of objects in different halos, 
resulting in P{k,z) = Pih(fc,z) -|- P2\i{k,z). The num- 
ber of pairs of galaxies within an individual halo is re- 
lated to the variance of the halo occupation distribution, 
fT^(M, z) = {Nga\{Nga\ — 1)), while the number of pairs 
of galaxies in separate halos is simply the square of the 
mean halo occupation number, N{M, z) = (Ngai) (HON, 
hereafter). We mo del the HON usin g a central-satellite 
formalism (see e.g. lZheng et al.l[2005t ): this assumes that 
the first galaxy to be hosted by a halo lies at its center, 
while any remaining galaxies are classified as satellites 
and are distributed in proportion to the halo mass pro- 
file. Different HODs for central and satellite galaxies are 
then applied. For central galaxies, the mean HON, A'con, 
is described by a step function such that halos above a 
minimum mass threshold M^un contain a single central 
galaxy and halos below this threshold contain no galax- 
ies. For satellite galaxi es, a power-la w in mass describes 
their mean HON (e.g.. iZehavi et al.l [2005) 



iVsat(M) = 



M 
Ml 



^cut /t^vit 


x'^ ■ 

A-min 


log(A/n,in/Mo) 


a 


1 

2 
3 
4 


16.3 
13.6 
11.5 
9.7 


11 Kn+OAO 

ii.OU_Q Q5 

11 rrn + OAO 

ii.OU_Q Q5 


'-'■^3-0.95 

1 icr-t-0.05 
-'-•-'-■'-0.75 


TABLE 3 



Best-fit values obtained for Mj^in and a for different 

CHOICES of the radius r^ut- THE MINIMUM-X^ (wiTH 10 
DEGREES OF FREEDOM) ARE ALSO SHOWN. 

where Mi is the mass-scale at which a halo hosts exactly 
one satellite galaxy (in additio n to the centra l galaxy). 
Both semi-analytic models (e.g.. lBerlind et aLir200Q) and 
hydrodynamical simulations (e.g.. IBerlind et al.l 120031 ) 
show that the distribution of galaxies within a halo is 
close to Poisson in the high-occupancy regime, i.e., when 
^sat 3> 1, and (strongly) sub-Poissonian in the low- 
occupancy regime. In order to agree with these results, 
satellite galaxies are assumed to be Poisson distributed 
at fixed halo mass. The distinction between central and 
satellite galaxies then automatically accounts for the sub- 
Poisson ian behavior of th e HOD in the low-occupancy 
regime ()Zheng et al.ll2005[ ). 
The 1- and 2-halo power spectra are 



Pihik, 



JT-halo {M, z) [27Vcon (A'/)7Vsat (M)mdM {k,z\M) 

Nl,{M)ulM{k,z\M)]dM/nl,i{z), 



lik, 



nhalo(M,2)iVgal(M,2;) 



b{M,z) UUM {k,z\M)dM 



(11) 



(10) 



The meaning of the symbols here is as follows: Pom is the lin- 
ear power s pectrum of dark matter, derived using the recipes 
of lEisenstein fc Hul l|l998i) for the matter transfer function; 



BLAST: Correlations in the Cosmic Far-Infrared Background 



9 



Whaio is the halo- mass function fsee ISheth et al.ll2001f ): 6 is 
the linear bias parameter; mdm is the normalized dark matter 
halo density profile in Fourier space; and rigai is the mean 
number of galaxies per unit comoving volume at redshift 2:, 



nhalo{M, z) 



M 



M 
Ml 



dM. (12) 



The expression for the 1-halo term implicitly assumes 
that the distribution of galaxies traces that of the 
dark matte r, for which wc have adopted the profile of 
INavarro. Fre nk. & White (1997, NFW), with the same con- 
centration parameter as Bullock ct al. ( 2001i ). Since the NFW 
profile formally extends to infinity, it is necessary to artifi- 
cially truncate the distribution at some radius, rcut- Typi- 
cally, this is chosen to be the virial radius of the halo; how- 
ever, this may not necessarily be realistic. We address this 
by first adopting the assumption that rcut = ^vir, and then 
exploring the consequences of relaxing that requirement, so 
that galaxies are allowed to lie further out. On large scales, 
where clustering is predominantly linear, udm ~ 1, so that the 
2-halo power spectrum simplifies to P2h ~ bifi{z)Pum(k, z), 
where bca{z) is the efi'ective large-scale bias. 



(-) 



ni,^i,{M,z)N^,i{M)h{M,z) dM/n^,i{z). (13) 



Our model has two free parameters, Mmin and a, which 
we vary through < a < 2 and 10 < log(Mmin/MQ) < 
16, with steps of 0.05 in both logM and a. The best- 
fit values of the parameters are determined through a 
minimization technique by fitting the observed power 
spectrum at each of the three BLAST bands simultane- 
ously. 

Throughout we assume that both M^m^ and a remain 
constant in time, although in principle they are functions 
of redshift. Whether these parameters evolve with red- 
shift would be difficult to constrain from our data alone; 
nevertheless, our assumption is consistent with what is 
observed fo r other class es of high- redshift sources (e.g., 
quasars, see lPorciani et al. 2004). For each M^i^-a pair, 
the mass-scale Mi is fixed by requiring that at every red- 
shift, z, the number density of the background sources 
derived from the halo model formalism matches that pre- 
dicted by the adopted source-population model, i.e.. 



Jo 



dN 
dS dz 



{S,z) dS = n^^i{z) dV^(z). 



(14) 



9. HALO MODEL FITS 



Under the assumption that the dark matter halos are 
described by an NFW profile truncated at the virial ra- 
dius, i.e., where rcut = ?'vir, and that galaxies within 
the halo trace the underlying dark matter distribution, 
the best fit to the observed angular power spectrum gives 
log(M„in/MQ) = 11.5^0:4 and a < 1.0, with x^n = 16-3 
with 10 degrees of freedom; i.e., the model is marginally 
consistent at the 2-a level. 

While this model is not formally ruled out by the 
data, as shown by a dotted line in Figure [6l it poorly 
reproduces the shape of the observed power spectrum, 
which exhibits a steeper slope. If we were to weight 
this model to fit the large scale power spectrum (e.g., 
k ^ 0.08 arcmin"^), then it would over-predict the power 



on small scales. Arguing that perhaps the model de- 
scribes the small scale, 1-halo term correctly, but is 
under-predicting the large-scales, is not a good descrip- 
tion because: (i) the 2-halo term is less sensitive to 
the underlying assump tions than th e 1-halo term (for 
a discussion see Tinker et al.l l2009f ) ; and (ii) the ob- 
served small-scale power spectrum is still too steep to 
be accounted for by the shallower 1-halo term. An- 
other possibility is that the redshift distribution of the 
background cumulative flux predicted by the adopted 
source count model is incorrect. In order to reproduce 
the observed shape of the angular power spectrum, the 
bulk of the background would have to o riginate from 
source s at z < 1, wh i ch is ruled out bv iDevlin et al.l 
(|2009f) . lMarsden et all (|2009[ ). and lPascale et al.l (|2009r 
However, a full investigation of the leeway in changing 
the redshift distribution (and degeneracies with other 
changes to the model) are beyond the scope of the present 
study. 

In light of this, we explore the possibility that the dis- 
crepancy between the predicted power spectrum and the 
observed one is related to the modeling of the 1-halo 
term. There are two obvious ways to modify the shape 
and normalization of the 1-halo power spectrum. One is 
to allow the dark matter halos, although still following 
an NFW profile, to be truncated at a scale rcut > ''i.ir- 
Thus, satellite galaxies are distribut ed over a larger vol- 
ume. This idea is not new; Magl iocchetti fc Porcianil 
(2003*) find that in order to adequately fit the 1- and 
2- halo term to the 2dF Galaxy Redshift Survey data-set 
it is necessary that the galaxies are allowed to reside out 
to 2 times the v irial radius. Further more, from semi- 
analytic models iDiaferio et al.l ()1999[ ) show that blue 
(and hence star-forming) galaxies tend to reside in the 
outskirts of their host halos, while red galaxies are found 
closer to the halo center. 

The second possibility is that the distribution of galax- 
ies within the halos does not follow that of the underly- 
ing dark matter. For example, a power-law distribution 
p(r) oc r"''' with 7 < 2 would make the 1-halo angu- 
lar power spectrum steeper than that predicted by an 
NFW profile. Simi lar arguments have been made by 
iWatson et al.l (|2009D . who found that by allowing the 
inner slope of the density profile to vary, they fit the 
small-scale clustering of luminous red galaxies quite well. 

Here we only explore the first possibility, examining 
'"cut — 1, 2, 3, and 4 x rvir. The results are shown in Fig- 
ure[6l and summarized in Table[3l The best-fit model an- 
gular power spectrum for the case rcut = 3 x Tvir is shown 
in Figure[7l while the corresponding HON and large-scale 
effective bias are shown in Figure [8] Note that while the 
value of the reduced Xmiu approaches unity for increasing 
values of rcut, the best-nt values of Mmin are negligibly 
affected by the changes, while a only marginally increases 
with increasing rcut- The effective mass of the halo, M^s, 
is the weighted mean over the halo-mass distribution: 

M,tf{z)^ [ nhaio(M,z)Ar(M,z)MdAf/ngai(z). (15) 

The results for the large-scale efi'ective bias b, and mass 
log(M/M0), at the respective medians of the redshift distri- 
butions of the sources contributing to the background in each 
band (see Figure [5} are 2.2 ±0.2, 2.4 ±0.2, and 2.6 ±0.2, and 
12.9 ±0.3, 12.8 ±0.2, and 12.7 ±0.2, at 250, 350 and 500 ^m, 
respectively. They are minimally affected by the change in 
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Fig. 7. — Power spectrum of background correlations from clustering of extra-galactic sources measured in the BLAST maps (circles with 
error bars) overlaid with the best-fit halo model (thick solid line), under the assumption that dark matter halos are NFW spheres truncated 
to 3x the virial radius, and the galaxies within the halo follow the underlying dark matter distribution. The shaded region shows the 99% 
confidence region in the Afmin-a parameter space. The dashed and the dot-dashed curves show the 1- and 2-halo contributions to the 
power spectrum, respectively. For comparison, the power spectrum obtained under the assumption that galaxies are unbiased tracers of 
the underlying dark matter distribution, i.e., Psoik, z) = PuM{k, z), is plotted as a lighter solid line. According to the model, the BLAST 
data occupy a range of angular scales which should be sensitive to both the linear and non-linear clustering terms. 
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U (M3„„) redshift z 

Fig. 8. — Left-hand panel: Halo Occupation Number (HON) as a 
function of the mass of the halo at a representative redshift z ~ 1.2, 
assuming the best-fit values Mniin = lO^"*^'^ Mq and a = 0.95. 
Right-hand panel: Redshift dependence of the large-scale effective 
bias, resulting from the best fit to the observed power spectrum 
of correlations due to clustering. The redshift interval over which 
68% of the signal originates, 0.7 ^ 2.5, is shown as a thick solid 
curve. 

fcut, changing by less than 8% for the bias, and 3% for the 
mass, over the full range of rcut expored. 

Figure shows the contribution to the total clustering 
power spectrum from sources in different redshift slices, 
within the assumed source population model. As expected, 
the power is dominated by the contribution from sources in 
the range 0.7 < z < 1.5, with an increasing contribution from 
sources at 2 = 1.5 — 3.0 for increasing wavelengths, as is ex- 
pected from Fig ure [5l T his is consi stent with the find ings of 
iDevlin et all (j2009); Ma rsden et al.l (2009); and Pasc ale et al.1 
([2009), who through stacking show that of the sources mak- 
ing up the CIB, the fraction at 2 > 1.2 increases from 40% at 
250 ^m, to 50% and 60% at 350, and 500 fitn. 

Finally, we can use the model to interpret the clustering 
in terms of a 3D spatial correlation length, ro. To do that, 
we Fourier transform the best-fit power spectrum to find 
the spatial correlation function, ^(r), from which ro is then 
simply the linear comoving scale at which the correlation 
function equals 1 at each redshift. This is a model-dependent 
approach to estimating ro, and as such should be considered 
an approximation. The more typical approach, which 



involves finding the angular correlation length, the redshift 
distribution, and deprojecting the signal by inverting the 
Limber equation, would result in very large uncertainties. 
The model-dependent ro is only mildly sensitive to the 
choice of rcut, varying by 10% over the full range. The 
model-dependent values for ro are illustrated in Figure [10] as 
a solid line with a shaded area representing 3-a uncertainties. 
The three BLAST points are plotted at the fiux-weighted 
median redshifts of the unique redshift distributions probed 
by the three bands (i.e., the distributions shown in Figure|4]), 
and should not be interpreted as the locations of all the 
sources contributing power in those bands. At these effective 
redshifts, (i.e., 2 = 1.60, 1.86 and 2.15), we find ro = 4.9, 5.0, 
and 5.2 ±0.6. h^^ Mpc at 250, 350 and 500 /im, respectively. 

10. DISCUSSION 
10.1. Comparison with other observations 

Comparisons with other measurements of clustering must 
be made and interpreted with care because not everyone uses 
the same definition of bias, the same parameterization for 
the halo model, etc. Nevertheless it is interesting to put our 
measurements into the context of the large body of literature 
on the clustering of galaxies selected in difi'erent ways, in order 
to understand how they might be related. 

In Figure \TU\ we compare the correlation lengths vs. red- 
shift of star-forming populations selected using a wide variety 
of techniques. It should be noted that in som e cases these 
techn iques select overlapping populations (see IReddv et al.l 
[2005I . for a nice discussion on the overlap between color se- 
lected samples). This list is by no means exhaustive, and 
is meant only to be an illustration. The reported values of 
ro were converted assuming a fixed slope 7 = 1.8, such that 
ro,i.8 = (ro,7) '^^'^ . The BLAST best-fit model estimates are 
shown as a line and shaded 3-cr confidence region. The three 
BLAST points are located at the median locations of their 
respective redshift distributions (see Figure In addition, 
the simulated clustering lengths of dark matter halos of given 
mass and redshift are measured from the Millennium Sim- 
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Fig. 9. — Contributions to the total clustering power spectrum from sources in increasing redsliifts slices. BLAST measurements and 
the best fit to the halo model (with rcut = 3 X rvir) are shown as circles with error bars and a solid line, respectively. Overlaid are 
the contributions from: dotted line (green), < 2: < 0.7; dashed line (red), 0.7 < 2 < 1.5; dot-dashed line (blue), 1.5 < z < 3.0; and 
triple-dot-dashed line (magenta), 2 > 3. It is clear that at 250 /jm, the bulk of the signal comes from galaxies in the redshift range 0.7—1.5, 
and that the contribution from galaxies in the redshift range 1.5—3.0 increases with increasing wavelength. 



ulation^ (|Springel et al.l [20051 ). by fitting a single power- law 
witli slope of -1.8 to the correlation function, and are shown 
as dotted lines. 

It is immediately clear that the galaxies which make up 
the background are not as strongly clustered as the more lu- 
minous sets of resolved sources (with the exception of those 
identified by their Lyman break, i.e., BM, BX and LBG). 
Furthermore, the strength of the clustering increases with 
incre asing luminosity (e.g., iGilli et al] 120071 : iBrodwin et al.l 
[2OOI 1. Since each of the techniques used to select the pop- 
ulations of galaxies that lie above the BLAST lines has an 
IR component, it is tempting to conclude that all of these 
populations contribute to the total submillimeter GIB. The 
relative contribution of each of th ese popula tions could be ex- 
plored through stacking, as Mars den et al.l ([2009) have done 
for BzKs. They found that although the BzKs make up about 
a quarter of the sources which completely resolve the GIB, 
they contribute ~ 32%, ~ 34%, and ~ 42%, at 250, 350, and 
500 /xm, to the total BLAST intensity. This indication that 
the resolved sources pick out parts of the background high- 
lights the complementary nature of the clustering measure- 
ments of the GIB and resolved sources in forming a complete 
picture of the environments of star-forming galaxies. 

Although the correlation length of the background galax- 
ies appears to change very little with redshift, the bias is a 
strong function of redshift (see Figure|8]). While both the bias 
and the correlation length are indicators of galaxy clustering, 
their behavior is not in contradiction because the clustering 
strength of the host dark matter halos is rapidly increasing 
with decreasing redshift as well, thus in this sense, it is the 
bias that is a more telling description of how star formation re- 
lates to structur e formation. This str ong evolution of the bias 
is confirmed bv lLagache et all (|2007l ). who found a redshift- 
independent bias parameter, b ~ 2.4, for the sources which 
make up the GIB at 160 ^m^, at 2 ~ 1. Our earlier result (see 

^ Note that Millenium Simulation cosmology is 
erg = 0.9 and Cm = 0.25, which has a minimal im- 
pact on the prediction. Catalogues can be found at 
http: / /www. mpa-garching.mpg.de/Millennium 

° Note that the value 6 = 1.7 reported in Lag ache et al.l {2007) is 
for (Tg = 1.1, and not eg = 0.8, as quoted in the text. The correct 
value of the measured bias parameter is 2.4 ±0.2 (Lagache, private 



SectionO, for galaxies which lie at higher redshifts, was fe ~ 4. 
Thus, in the scenario of a strongly evolving bias parameter, 
from anti-biased in the local Universe, to highly bia sed at 
z ~ 1 and beyond (e.g.. lSheth eraL|[200^: lElbaz et al.ll2007, ). 
our result is consistent. 

10.2. Clustering in Context 

ILe Floc'h et al.l (|2005l ) show that IR-luminous galaxies 
(LIRGs and ULIRGS) represent ~ 70% of the IR energy den- 
sity, and are responsible for most of the sta r formation, at 
z ~ .5-1.0 and beyond. Stacking work (e.g.. iMarsden et al.l 
I2009D shows that most of the objects responsible for the GIB 
are fainter than the flux density limit of the BLAST catalogs 
(j Pevlin et al.|f2009( ), corresponding to LIRG-like luminosities 
for those sources. Therefore, by identifying the locations of 
active star formation, we are also identifying the locations 
of the formation of the majority of stars in the present day 
Universe. With this in mind we ask: where are active star- 
forming galaxies preferentially found through cosmic time? 
From Figure 1101 it appears that the most active star-forming 
galaxies are found occupying halos whose mass increases from 
roug hly 10^2 to lO^^ '^ Mq over the redshift range z = to 
2, after which it appears to remain roug hly constant. Thi s 
appears to be consistent with downsizing (jGowie "eFallligQel ). 
the scenario in which the sites of most active star formation 
shifts to ever larger galaxies at higher redshifts. However, 
according to the model, which is only constrained for red- 
shifts greater than ^ 1.1, the galaxies which make up the 
background do not exhibit the same trend. While the bias 
is still a strong increasing function of redshift, the clustering 
strength of the host halos remains roughly constant, corre- 
sponding to typical host halos that become slighly smaller. 
Since this is very much a model-dependent claim, it may be 
indicitave of a flaw in the model (perhaps assuming Mmin-ct 
remain constant with z is incorrect); future studies should 
clarify this picture. 

A striking feature is the sharp cut-off at M > 10^"^'^ Mq, 
which appears to hold out to z ~ 2.5. As was pointed 
out by IBrodwin et al.l l|2008D . this appears to be inconsis- 
tent with models which claim that star formation should be 

communication) . 
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Fig. 1 0. — Comoving c orrela tion leng th vs. redshift for star-forming galaxies selected with a variety of techniqu es. Ot her data taken from: 
IRAS -'Saunders et al.' fT99l); SMG -IWebb et aP (f2003); 'Blain et alJ (!2Q0l); B2 and B3 -IFarrah et"an (f20M); IR -'Magliocchetti et alj 
(12007, 2008); Gilli ct al. (2007); DRG - Grazian ct al. (2006); Quadri et al. (2008) ; DOG - [Brodwin ct al. (2008); BzK - Blanc et aQ 
j2008); Hartlcv ct al. (2008); BM and BX - Adclbcrgcr ct aj^ 12005 ); LBG - Giavalisco et al.l ( [1998 ); Adelberger et al., t2005); and UVLG 
— [Basu-Zvch ct al. ( 2009) . The dotted fixed mass lines show the predicted clustering lengths of halos of a given mass at a given redshift, 
found with the Millennium Simulation I jSpr ingcl ct al. 20011). The model-dependent BLAST values for ro are shown as a solid line with a 
shaded area representing 3-a- uncertainties. The three BLAST points are plotted at the median redshifts of the distributions from which 
the signal originates (i.e., the distributions shown in Figure [Sj. The ranges from which 90% of the power originates are il lustrated as 
corresponding colored lines. In the context of the model, the clustering strength of BLAST galaxies is compatible with the IWebb et al.l 
(|2p03) and Blain ct al.. (,2004) estimates for clustering of submillimeter galaxies, but less strong than the that of other resolved populations 
of galaxies. 



quenched in hal os with masses greater t han 10 due 
to shock heating (Birnboim & Dckcr2003VDekcI fc Birnboiml 
|2006; Cattanco ct al. 2008). Dckcl ct al. (2009) attempt to 
resolve this dilemma with a model where cold streams pene- 
trate the shock-heated media. On the othe r hand, if the shock 
radiu s roughly follows the virial radius (|Birnboim fc Dekell 
I2003D . then finding satellites actively forming stars outside of 
the shock-heated volume would satisfy both the model and 
the observations. 

11. CONCLUSIONS 

We report the detection of correlations in excess of Poisson 
noise in the GIB, over scales of approximately 5-25 arcmin, 
with BLAST at 250, 350, and 500 /im, at a level with respect 
to the GIB of AJ/I = 15.1 ± 1.7%. 

The GIB is made almost entirely out of individual sources 
distributed over a wide range of redshifts. We find that within 
the context of a reasonable model for the source popula- 
tion, the signal originates from galaxies in the redshift ranges 

1.3 < z < 2.2, 1.5 < z < 2.7, and L7 < z < 3.2, with median 
redshifts z = 1.61, 1.88 and 2.42, at 250, 350 and 500 ^m, 
respectively. Fitting to the linear theory power spectrum, we 
find that the BLAST galaxies responsible for the GIB fiuc- 
tuations have a bias parameter, 6 = 3.8 ± 0.6, 3.9 ± 0.6 and 

4.4 ±0.7. 

We further interpret our results in terms of the halo model. 
We find that the simplest prescription does not fit very well. 
One way to improve the fit is to increase the radius at which 
we artificially truncate dark matter halos to well outside the 
virial radius. This may imply that the star-forming galax- 
ies that we are seeing at z ~ 1 are preferentially found in 
the outskirts of groups and clusters. This is consistent with 
related phenomena that have been observed at other wave- 



lengths (|Magliocchetti et al.1 120041: IMarcillac et al.] 120071 ). as 
well as in simulations (jPiaferio et al.lll999l ). 

For a HOD with 'satellite' galaxies occupying halos out 
as far as rcut = 3r-vir, we find parameters log(Mniin/MQ = 
11.5*^0.1, and a = l.lj^oi' resulting in efi'ective biases 
6efi = 2.2 ± 0.2, 2.4 ± 0.2, and 2.6 ± 0.2, and effective masses 
log(Mcff/MQ) = 12.9 ±0.3, 12.8 ±0.2, and 12.7 ± 0.2 at 250, 
350 and 500 /xm, corresponding to spatial correlation lengths 
of ro = 4.9, 5.0, and 5.2 ± 0.7 h"^ Mpc, respectively. 

In the context of the model, we see that star formation is 
highly biased at z ^ 1, unlike in the local Universe, where 

analogous galaxy populations, such as IRAS galaxies, are 
found to be mildly anti-biased (e.g., M^s 10^ M© and 

h = 0.86, ISaunders et"aIlfT993 ). 

We find relatively small values for Oq, further confirming 
that the sources which make up the GIB are distributed over 
a wide range of redshifts, which has implications for planning 
future submillimeter clustering measurements. For example, 
as Knox et al. (2001) have argued, to match the precision of 
the bias measurement made from correlations in the back- 
ground by using discrete sources will be very challenging. To 
achieve ~ 5% accuracy would effectively require thousands 
of sources with exact redshifts, over tens of square degrees. 
When redshifts are only approximately known, this increases 
to hundreds of thousands of sources over hundreds of square 
degrees, which is only possible with instruments whose beams 
are smaller than 5 arcsec. Thus, while measuring the cluster- 
ing from resolved sources has numerous advantages - for ex- 
ample studying the clustering properties of a subset of galax- 
ies - accurately measuring the large scale bias is best achieved 
through correlation analysis of the background fluctuations. 
This will be a focus of future studies with BLAST, as well as 
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with the Herschel and Planck satclhtcs, and SCUBA-2. 
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